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Abstract 

In this paper, we apply the recently developed ab initio renormalized excitonic method (REM) to 
the excitation energy calculations of various molecular aggregates, through the extension of REM 
to the time-dependent density functional theory (TDDFT). Tested molecular aggregate systems 
include one-dimensional hydrogen-bonded water chains, ring crystals with tt-tt stacking or van-der 
Waals interactions and the general aqueous systems with polar and non-polar solutes. The basis 
set factor as well as the effect of the exchange-correlation functionals are also investigated. The 
results indicate that the REM-TDDFT method with suitable basis set and exchange-correlation 
functionals can give good descriptions of excitation energies and excitation area for lowest electronic 
excitations in the molecular aggregate systems with economic computational costs. It's shown that 
the deviations of REM-TDDFT excitation energies from those by standard TDDFT are much less 
than 0.1 eV and the computational time can be reduced by one order. 

Keywords: fragment-based quantum chemical methods; time-dependent density functional the- 
ory (TDDFT); supra-molecule; excitation energy; low-scaling; Frenkel exciton model 
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I. INTRODUCTION 



Molecular aggregates are coupled clusters of small molecules with intermolecular separa- 
tions typically close to individual molecule size, for example, the biological photosynthetic 
light harvesting system, the organic semiconductor crystal or the solute dissolved in sol- 
vents. Moreover, the ability of converting solar light into electrical or chemical energy in 
these systems through photosynthesis or photoelectric conversions motivates the study of 
the electronic excited states of molecular aggregates. However, the theoretical characteri- 
zation of these properties is often challenging to unravel due to their relatively large scales 
and complicated environments. Among the current popular quantum chemical methods for 
calculating electronic excited states, time- dependent density functional theory (TDDFT) 
is mostly widely used due to the good balance between the accuracy and computational 
cost£~— . It is well known that standard approximate exchange-correlation functionals used 
in DFT or TDDFT will underestimate the excitation energies for Rydberg states and charge- 
transfer states as well as extended 7r-conjugated systems and weakly interacting molecular 
aggregates. Such drawbacks are due to the fact that those functionals do not exhibit the 
correct ^ asymptotic behavior and can not capture long-range correlation effects. Recent ef- 
forts have offered possibilities to account for long-range corrections and dispersion effects by 
the newly developed exchange-correlation functionals with long-range corrections 4 ^— and/or 
dispersion corrections^—. However, the applicable system size for excited state quantum 
chemistry calculation is still limited to a few hundred atoms at the most. Since the N 3 ~ 4 
scaling (N is the size of system) of the TDDFT—, the application of TDDFT to very large 
systems is still challenging. 

In order to reduce the computational scaling in TDDFT, many theoretical approaches^— 
based on the local correlation approximation have been suggested. Chen and co-workers^ 
developed a linear-scaling time- dependent density functional theory (TDDFT) algorithm 
using the localized density matrix (LDM) and in an orthogonal atomic orbital (OAO) 
representation. Yang and co-workers^ extended this formalism and suggested reformu- 
lating TDDFT based on the non-orthogonal localized molecular orbitals (NOLMOs)^ 1 ^. 
Casida and Wesolowski proposed the TDDFT within the frozen- density embedding (FDE) 
framework—, and Neugebauer and co-workers extended this approach with coupled electronic 
transtions 4 ^ 4 ^ and made applications of this approach to many interesting systems 4 ^— like 
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light-harvesting complexes in biomolecular assemblies^. Recently, based on the fragment 
LMOs that derived from capped fragments, Liu and co-workers^ suggested a new linear- 
scaling TDDFT method and successfully applied it to several large conjugated systems. 

Considering the weak interactions between the molecular units, using a "divide and con- 
quer" idea to treat the excited states of the aggregated systems may be a worthwhile 
attempt^. In the fragment molecular orbital (FMO) method proposed by Kitaura and 
Fedorov^ 1 ^, the whole system can be divided into small fragments, and the total proper- 
ties can be well estimated by the corresponding monomers, dimers, etc. Recently, using the 
FMO scheme, FMOx-TDDFT— ~— (x means n-body expansion) with analytic gradients have 
been developed and can give good descriptions for solvated and bio-chemical systems. Mata 
and Stoll^ also developed an improved incremental correlation approach for describing the 
excitation energies, with the inclusion of a dominant natural transition orbitals into selected 
excited fragment. However, these methods may lose efficiency when dealing with general 
systems which have multiple or uncertain excited regions. 

For general systems, the Frenkel exciton modePi may be used as an alternative subsys- 
tem strategy and this model has been first applied to molecular crystals and subsequently 
extended to aggregates^— . In its original form, the Frenkel exciton Hamiltonian describes 
a weakly interacting ensemble of two-level systems by 



where indices % and j label "blocks" (molecules), bf (hi) are the creation (annihilation) 
operators of an excitation on block i, fij is the excited-state transition energy of block i, and 
Sij is the interaction, or coupling, between blocks i and j. The direct product of eigenstates 
of isolated blocks forms a convenient basis set for the global excited states. However, the non- 
diagonal couplings (sjj) between different blocks, if calculated within the truncated Hilbert 
space directly, involve only the electrostatic interactions, and the description of excited 
states with such an approximation will fail for the systems in which the quantum dispersions 
dominate the inter-molecular interactions. Improvements have also been proposed through 
empirical corrections or exchange-correlation potentials^ - — . 

In the recent years, the contractor renormalization group (CORE) method^ 1 ^ and also 
the real-space renormalization group with effective interactions (RSRG-EI)^ provided a 
novel kind of subsystem methods for describing excited states of large systems, in which 
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the general excited state is assumed to be an assembly of various block excitations (both 
single block excitations and multiple block excitations), and the interactions between ad- 
joined blocks are taken into account through Bloch's effective Hamiltonian theory^ 1 ^. The 
basis set used in CORE or RSRG-EI is similar to that of Frenkel exciton model, but the 
non-diagonal coupling terms in CORE or RSRG-EI include not only the electrostatic inter- 
actions but also the quantum exchange contributions through the super-block calculation 
and the followed projection of the super-block wave-function onto the blocks' direct product 
basis, in contrast with those in traditional Frenkel exciton models which involve only the 
electrostatic interactions. In 2005, Malrieu and co-workers^ made the further simplification, 
approximating the excited states of the whole system as only the linear combination of var- 
ious single block excitations, and proposed it as the renormalized excitonic method (REM). 
Recently, we and our co-workers^^ applied the REM into ah initio quantum chemistry and 
successfully combined it with various ah initio methods like full configuration interaction 
(FCI), configuration interaction singlet (CIS), and symmetry adapted cluster configuration 
interaction (SAC-CI). Good descriptions for excited states and ionized states of hydrogen 
chains and polyenes as well as polysilenes with economic computational costs have been 
achieved with both orthogonal localized molecular orbitals (OLMOs) and block canonical 
molecular orbitals (BCMOs)^. 

In this paper, we combine REM with TDDFT and extend the REM calculated systems 
from the linear molecule to various molecular aggregates. Testing systems include hydrogen- 
bonded H2O molecular chains, ring crystals like vdW interacted H2O rings and n-n stacked 
C2H4 rings, 2-D benzene aggregates, as well as aqueous systems with polar and non-polar 
solutes. The REM-TDDFT wavefunction as well as the basis set and functional factors are 
also discussed. The structure of this paper is designed as: In Sec. II, technical details of 
the method are introduced; in Sec III, we present calculated results of various molecular 
aggregates using REM-TDDFT and make comparisons with standard TDDFT calculations; 
and finally, we summarize and conclude our results in Sec. IV. 

II. METHOD 

The REM method is a type of fragment-based method. In REM, the whole system can 
be divided into many blocks (usually tens or hundreds of), as illustrated in Fig. [U Here the 



4 



I, J, K, L are the block-monomers just like in Frenkel exciton model, and additionally, the 
adjacent monomers form the dimers. 



dimer-JK 



monomer-i 



monomer-J 



dimer-IJ 



monomer-K 



monomer-L 



FIG. 1: The partition of the whole system into various blocks 

In our previous work—, we gave a detailed description of constructing REM Hamiltonian 
with BCMOs. Under BCMOs, the orthogonality only exists in the intra-blocks orbitals but 
not in the inter-blocks orbitals, and we used an approximate projector (Pq) to unravel the 
non-orthogonality situation^. Since this strategy is universal, we can use the subsystems' 
Kohn-Sham molecular orbitals (KS-MOs) instead of the BCMOs, and here we briefly in- 
troduce our REM-TDDFT strategy as below (up to 2-body interactions). For more details 
about the method, we refer the readers to our previous work 78 . 

1) Calculate each block- monomer by TDDFT to get KS-MOs, as well as the eigenstates 
for ground state (ip°) and an excited state (i/j*). In our REM-TDDFT strategy, we assume 
the if) state is closed-shell ground state, and ip* state is constructed by the excitation 
components in TDDFT. Then the basis functions (\^rem)) of the model space are formed 
by 

N N 

i*™>=Ei^>=E[i^>ni^>] ( 2 ) 
i=i i=i 

where N is the total number of block monomers and the /, J means the i-th, J-th block- 
monomer, respectively. The corresponding projector (Pq) by 

N N 
Po = \^REm)S- 1 (^ REM \ = £ miS' 1 )!!' 

1=1 /,=1 (3) 

N N W 

= £[i^> n i$>](3r»v ek#i n <^°'0 
i=i jjki r=i j'jkr 

where S m is the overlap matrix between model space basis functions. 

Here we should mention that the antisymmetric property is satisfied in this product 
basis (Eq. |2|). Although the MOs (or the electron density) of each monomer is totally 



localized on its own, quantum mechanics (QM) calculations of the dimers (or trimer, etc.) 
can redistribute the electron density^, which is crucial to discribe the charge transfer type 
excitation. Owing to the multi-mer's contribution, REM can give reasonable descriptions 
for various types of the low-lying excitations. 

2) Obtain two lowest excited states of dimer (ipjj and ijjfj) and use these eigenstates to 
form target states 

i*?j> = wti n io + iv#> n \^k) (4) 

Following with projecting the target states onto model space via projector (P ) 



N N 

7iTr* llTf *^ 



(5) 

N N K ' 



r=i i=i 

Here we can denote the E/Li^m 1 ) n' (^*r\^*u)] as matrix C . How to solve C may be 
the most complicated part in the ab initio REM strategy, one can refer to the Ref.— for 
detailed illustration. Nevertheless, we should also mention that two excited states are chosen 
here, as a result of only one excited state is kept in each monomer; if one more excited state 
were kept in monomer, then four excited states should be appropriate. 

3) Once we get the Co matrix, usually the orthogonalization process should be applied to 
Co to get a new set of coefficients C, in order to obtain the Hermitian Hamiltonian. Then 
the expression of dimer-/J' effective Hamiltonian can be written in the matrix form 

Hf/ = (C + )- 1 e /J C- 1 (6) 

where the eu are the two lowest excited states energies of dimer-JJ. And the interactions 
between monomer-I and monomer- J can be acquired by 

H e /J = Hf/ - (Hf? + Hf) (7) 

Nevertheless, we should mentioned that the dimensions of various effective Hamiltonians 
are the same, and they are all equivalent to the number of REM basis. Here we take the 
H e /J (Eq. EJ) as example, as shown in Fig.[2j It could be found that the construction of H e /J 
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is under the whole REM bases while only two lowest excited states energies e*fj and e}j are 



used, then the HjJ will be a N x N matrix. When turn to (H^ 



eff 



Hf),the (eJ + 4) 



and (e5 + are used instead of efj and Sjj, and the dimension is also N x N. 



H eff = 
11 u 



* 



{vyy K -\ 




FIG. 2: Schematic illustration of the construction of H 



eff 
I J 



4) After the determination of the effective Hamiltonians for various dimers, the renormal- 
ized Hamiltonian for the whole system can be obtained according to the following expressions 



N 



reff 



i=i i>j 

and finally can be solved as the generalized eigenvalue problem, 

R eff C eff = S C eff E 



(8) 



(9) 



where the eigenvalues E are the excited state energies, and eigenstates C e ^ corresponds 
to contributions that the excitations occur in every blocks. In order to get the excitation 
energies, additional step should be applied to subtract the ground state energy with 2-body 
expansion (E — E ), where 



N 



i i>j 

Since the dimension of the H e ^ is equivalent to the number of REM basis, which is usually 
less than one thousand, the Jacobi Method can be used to diagonalize the H e ^ of the whole 
system. 



III. COMPUTATIONAL DETAILS 



We use our own code to implement the REM strategy. The preliminary TDDFT cal- 
culations on various blocks are implemented by GAUSSIAN— or GAMESS22.. In the 1-D 



H2O chains systems, we use GAUSSIAN09 to do the subsystem TDDFT calculations on 
the monomers and dimers. In this test, both of the adjacent and the separate dimers are 
considered. In the other systems, if no extra illustration, the modified FMO subroutine in 
the GAMESS package is used to automatically select dimers and do the subsystem calcula- 
tions. In our calculations the electrostatic potential (ESP) terms^&Si^ and other correc- 
tion terms are currently not considered, meaning only the original TDDFT calculations are 
implemented. 

IV. RESULTS AND DISCUSSION 
A. 1-D H 2 chains 

The model 1-dimensional water molecular aggregates are chosen as the starting test sys- 
tems. The geometrical configurations are represented in Fig. [3J There the typical hydrogen 
bond length (1.85A)^ and other two spacings (1.50A and 2.20A) are chosen. The O-H 
bond length is fixed as 0.9584A, and the angle of H-O-H is fixed as 104.45°. This starting 
system is simple and clear, and it should be a ideal model when we implement our starting 
REM-TDDFT calculation and do detailed analyses. 

First, long-range corrected exchange-correlation functional LC-BLYP and Pople's 6- 
31+G* basis functions are used here to perform the REM-TDDFT calculations. And in 
order to estimate the accuracy of REM-TDDFT, the standard TDDFT calculations are 
performed by GAUSSIAN09^. When performing the REM-TDDFT, two different fragmen- 
tation schemes are used here: fragmentation-A with one water molecule as one monomer, two 
water molecules as one dimer; fragment at ion- B with two water molecules as one monomer, 
and four water molecules as one dimer. In fragmentation-A, each monomer keeps the ground 
state and one excited state. The excited state can be Si or 7\, depends on which state you 
want to calculate (of the whole system). The dimers here keep the lowest Si and S% (or 
Ti and T 2 ) states. In fragmentation-B, the monomers and dimers contain double number 
of water molecules comparing to those in fragmentation-A, then there monomers keep one 
more excited state (S 2 or T 2 ), and dimers keep two more excited states (S3, S4 or T 3 , T 4 ). All 
the electronic structure calculations on various monomers and dimers are also implemented 
by GAUSSIAN09. 

The REM-TDDFT results are summarized in Table. [H Let's start from the Si state and 
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r t jX r r r r i 

Spacing : R=1.50, 1.85, 2,20 (A) 

FIG. 3: The 1-D water molecular aggregates 

take the 1.85 A spacing case (the typical hydrogen bond length) as an example. In REM- 
TDDFT with fragmentation-A, the standalone H2O monomer's excitation energy for S± is 
8.098 eV, and the standalone (H20)2 dimer's excitation energy for S± is 7.887 eV. Here the 
results of REM-TDDFT in (H 2 0) 8 , (H 2 0) 16 and (H 2 0) 24 are all 7.887 eV for Si states. 
These values are agree with the full TDDFT quite well, with the error of about 0.01 eV. In 
REM-TDDFT with fragmentation-B, the standalone (H 2 0) 2 monomer's excitation energy 
for Si is 7.887 eV, the (H 2 0) 4 dimer's Si is 7.895 eV. There the corresponding REM results 
are all 7.895 eV in the different size of water aggregates, and match quite well with the 
full TDDFT values (7.899 eV). The Si states in 1.50A and 2.20A spacing cases have the 
behaviours similar with the 1.85A case, and the errors are all less than 0.04 eV. One may be 



TABLE I: Calculated singlet /triplet excitation energies (in eV) by standard TDDFT with 
LC-BLYP/6-31+G* and related excitation energy differences between REM-TDDFT and 
them 



System 




Si 






s 2 






Ti 






T 2 




TDDFT 


REM' 4 


REM S 


TDDFT 


REM A 


REM S 


TDDFT 


REM A 


REM 3 


TDDFT 


REM- 4 


REM S 


2.20A 


























(H 2 0) 8 


7.942 


+0.018 


+0.024 


7.958 


+0.017 


+0.024 


7.198 


+0.021 


+0.018 


7.212 


+0.033 


+0.013 


(HaO)ia 


7.930 


+0.030 


+0.036 


7.936 


+0.033 


+0.046 


7.188 


+0.031 


+0.028 


7.192 


+0.050 


+0.027 


(H 2 0)24 


7.928 


+0.032 


+0.038 


7.931 


+0.038 


+0.051 


7.185 


+0.034 


+0.031 


7.188 


+0.053 


+0.030 


1.85A 


























(H 2 0) 8 


7.899 


-0.012 


-0.004 


8.000 


+0.092 


+0.073 


7.192 


0.000 


0.000 


7.284 


+0.129 


+0.033 


(H 2 0)i 6 


7.899 


-0.012 


-0.004 


7.978 


+0.111 


+0.094 


7.191 


+0.001 


+0.001 


7.263 


+0.148 


+0.049 


(H 2 0) 24 


7.899 


-0.012 


-0.004 


7.973 


+0.116 


+0.099 


7.191 


+0.001 


+0.001 


7.258 


+0.153 


+0.053 


1.50A 


























(H 2 0) 8 


7.575 


-0.004 


-0.016 


7.924 


+0.048 


+0.020 


6.955 


+0.029 


-0.002 


7.288 


+0.169 


+0.035 


(H 2 0)ie 


7.576 


-0.005 


-0.017 


7.890 


+0.079 


+0.035 


6.955 


+0.029 


-0.002 


7.252 


+0.196 


+0.057 


(H 2 0) 24 


7.577 


-0.006 


-0.018 


7.883 


+0.083 


+0.039 


6.955 


+0.029 


-0.002 


7.245 


+0.201 


+0.062 



A REM-TDDFT with H 2 as one monomer, (H 2 0) 2 as one dimer 
B REM-TDDFT with (H 2 0) 2 as one monomer, (H 2 0) 4 as one dimer 
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confused about why the S\ values in REM-TDDFT are not changed with the elongation of 
water chain. Here we take (H20)s for example to explain. From the result of full TDDFT 
calculation, some important orbitals are shown in Fig. H] and some important configurations 
(excitation components) of Si are listed in Table. [Ill It could be found that, the excitation is 
mainly from the HOMO (40-th) orbital to various unoccupied orbitals. Combining with the 
Fig. HI it could be found that the electron mainly from left-most H20's P z orbital, excites 
to various unoccupied orbitals (no P z components) for Si state. Since the excitation mainly 
stems from the left-most water, the elongation of water chain will not change the Si excited 
energy much. We also show the Si wave functions obtained from the REM-TDDFT in the 
left part of Fig. We could find from the coefficient analysis of REM wave function: the 
Si excitation is mainly contributed by the excitation from the left-most water, and a little 
component from the second left water. This picture is matching well with the full TDDFT 
calculation. Since there are only little components from the water monomers which are not 
belonging to the "left two" , the excitation energy for the Si states will only slightly change 
with the increasing chain length. 

J 9 J 9- 9 J 9 J 9 J f < 

40 HOMO 




?- p • * * r< 44 *JMtf *4 * 
* #99 * t* f t* 45 ^sf 

* - >■ >■ - 



FIG. 4: Some important frontier orbitals in (H20)g chain 



TABLE II: Some important coefficients in the 5*1 TDDFT amplitudes of (H20)s 

Excitation mode Amplitudes Excitation mode Amplitudes Excitation mode Amplitudes 
40 -> 45 -0.31473 40 -> 44 0.30290 40 46 -0.25920 

40 -> 43 -0.19455 40 -> 47 -0.19361 39 43 -0.12090 

39 -> 44 -0.10897 
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FIG. 5: Calculated contribution coefficients of various local excitations by REM-TDDFT for S\ 
and S2 

REM calculations for the triplet excited states usually give less error than those for the 
singlets, and the triplet excitations usually have a more localized character than the singlet 
excitations^. In Table. HI we can find the numerical accuracies for Ti states are as well as 
those of S\ states. In 1.85A spacing cases, the errors are only about 0.001 eV, no matter what 
fragmentation schemes we use. In 1.50A spacing, fragment at ion- A has errors of about 0.03 
eV. When enlarging the fragment units, the errors reduce to 0.002 eV. The errors increase 
to around 0.031 eV when the spacing turns to 2.20A (similar in Si case). This is due to the 
fact that the near-degeneracy problem for the very weak coupled systems will decrease the 
accuracy of REM calculations^. 

When turn to higher excited states, one could find the REM-TDDFT method can also 
give relative good descriptions. Here we take the S 2 state in 1.85 A separation case as an ex- 
ample: In REM-TDDFT with fragmentation-A, the deviations between REM-TDDFT and 
full TDDFT are 0.092 eV, 0.111 eV and 0.116 eV for (H 2 0) 8 , (H 2 0)i 6 , (H 2 0) 24 , respectively. 
These deviations are larger than those in Si states (about 0.01 eV) case. These deviations 
will turn to small when using larger monomers and dimers. With fragment at ion- B, they are 

TABLE III: Some important coefficients in the S 2 TDDFT amplitudes of (H 2 0) 8 
Excitation mode Amplitudes Excitation mode Amplitudes Excitation mode Amplitudes 



38 -> 42 


0.31013 


36 - 


» 42 


0.29387 


35 -> 43 


-0.17889 


36 41 


-0.16899 


35 - 


» 41 


-0.16083 


36 -> 44 


-0.15601 


38 -> 43 


0.15096 


35 - 


» 45 


-0.12300 


34 -> 41 


0.11775 


34 -> 44 


-0.11517 


34 - 


» 48 


0.11011 
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0.073 eV, 0.094 eV and 0.099 eV, but still larger than those in Si case. One can also find 
the same trends in other spacing cases or in T 2 excited states. Why are the errors in S2 
(T 2 ) larger than those in Si (71)? In general, the accuracy of higher states depends on the 
lower states during the diagonalization process, therefore the error in the latter promotes 
a larger error in the former—. There we can also use wave function analysis to illustrate 
this phenomenon in the REM-TDDFT calculation. In Table. IHfl we list some important 
(largest) excitation components in S2 of (H 2 0)s. It could be found that excitation in S2 can 
originate from many orbitals (34, 35, 36 and 38), means this excitation may be contributed 
by many H 2 monomers. In REM strategy, we only consider the monomers and dimers, 
this hierarchical structure will lessen the change of domains on excitation and also affect the 
accuracy^. In the previous Si state, the excitation mainly located on the edge, then the less- 
ened change of domains may not affect the result much. However, in S2 state, the excitation 
range from more H 2 units, then the lessened change of domains will affect the result more 
than in the Si state. Generally speaking, the higher excited states would own more excited 
regions, then the results will have larger errors. There we also show the REM-TDDFT S2 
wave function in the right part of Fig. |5j It could be found that our REM-TDDFT gives a 
normal distribution wave function with the peak in the middle of water chain. This picture 
agrees well with the full TDDFT calculation. 
B. Ring molecular crystals 

Now, let's turn to ring molecular crystals. Here we choose H 2 ring crystals dominated 
by vdW interacting and C 2 H4 ir-rc stacked ring crystals to test. These systems can be 
seen as simplified systems with periodic boundary conditions. In this part, we also check 
whether the error is affected by the basis set and the DFT functional. The geometries and 
the detailed descriptions of these two types of ring crystals are refer to Ref.— , water is put 
anti-parallel and ethylene parallel to each other. The alternating arrangement of waters 
would be favourable from the dipole-dipole interaction between the water monomers. The 
inter-water distance is 3. OA, near the 3.26A in which the ground state is attractive and has 
a minimum reflecting on the dipole-dipole interaction^. The inter-ethylene distance is 4.5A, 
near the 4.90A in which the inter-ethylene electron transferred state is attractive and has a 
minimum at around 4.90A. In this inter-ethylene distance, there are two types of excitations: 
one is ir — > 7r* excitations within each monomer; the other is electron-transfer type n — > n* 
excitations between monomers^. Parts of them are shown in Fig. [7J 
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R=3.0A 



R=4.5A 



• 



HO 10-mer 



C 2 H 4 50-mer 



FIG. 6: The geometries of the two types of ring crystals 



The long-range corrected functional (LC-BLYP) with three different basis functions (6- 
31G, 6-31+G* and 6-31 1++G**) are used here to perform the REM-TDDFT calculations 
and standard TDDFT. The subsystem TDDFT and standard TDDFT calculations are im- 
plemented by GAMESS^. The results are listed in Table. HVi In this table, we use two 
fragmentation schemes: the former (REM A ) is one H 2 (or C2H4) unit as one monomer, 
two H2O (or C2H4) units as one dimer; the latter (REM B ) is two H2O (C2H4) units as one 
monomer, then four H 2 (C2H4) units as one dimer. Each monomer keeps one ground state 
(So) and one excited state (Si), each dimer keeps two lowest excited states (Si, S 2 ). It could 
be found from the table that with former fragmentation scheme, the typical deviations in wa- 
ter ring systems between REM-TDDFT and standard TDDFT are about 0.06 eV in 6-31G, 
0.14 eV in 6-31+G* and 0.11 eV in 6-311++G**, respectively. When using latter fragmen- 
tation scheme, those derivations turn to -0.01 eV, -0.06 eV and -0.07 eV, corresponding. It 
could be found that, no matter what fragmentation scheme is chosen, the errors for the larger 
basis sets are only slightly larger than the smaller basis sets on the average. The similar 
tendency can also be found in the ethylene ring systems. In general, the subsystem methods 
have a somewhat larger errors with extensive basis sets for the excitation energy. This is 
because the interactions betweens fragments will be enforced in extensive basis sets, such as 
the exchange-repulsion and charge transfer, however, there the REM-TDDFT method can 
only recover the interactions from two body level, it isn't enough. 

Next, it is also of interest to see if the error is affected by the DFT functional. There we 
compare the Si excitation energies using BLYP, B3LYP with 6-31+G* basis, and also add 
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R=4.5A 
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C 2 H 4 50-mer 



FIG. 7: The geometries of the two types of ring crystals 

the LC-BLYP data in Table. IIVI There we also use the two fragmentation schemes as above. 
The results are summarized in Table. |V] We can find that the REM with long-range corrected 
functional LC-BLYP give a better description than the pure functional BLYP and the hybrid 

TABLE IV: Calculated S± excitation energies (in eV) by standard TDDFT and related 
excitation energy differences between REM-TDDFT results and them with different basis 
sets in the ring molecular crystals 



System 

H2O ring 
(H 2 O)i 
(H 2 O) 20 
(H 2 O) 50 



LC-BLYP/6-31G 



LC-BLYP/6-31+G* LC-BLYP/6-311++G 



TDDFT REM A REM B TDDFT REM A REM 5 TDDFT REAL 4 REM B 

6.763 +0.056 -0.003 6.840 +0.076 -0.067 5.857 +0.078 -0.069 

6.758 +0.066 -0.045 6.837 +0.153 -0.054 5.841 +0.113 -0.081 

6.756 +0.065 -0.005 6.845 +0.144 -0.065 5.836 +0.111 -0.084 



C 2 H4 rim; 

(C 2 H4)io 

(C 2 H4) 2 o 

(C 2 H4)5Q 



8.132 +0.007 -0.042 7.108 -0.022 -0.071 6.886 -0.010 -0.060 
8.119 +0.007 -0.043 7.080 -0.031 -0.075 6.843 +0.022 -0.024 
8.117 +0.005 -0.046 7.074 -0.034 -0.080 6.845 +0.018 -0.027 



A REM-TDDFT with H 2 or C 2 H 4 as one monomer 

B REM-TDDFT with (H 2 0) 2 or (C 2 H 4 ) 2 as one monomer 
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functional B3LYP. When using BLYP or B3LYP functional, the typical deviation in REM" 4 
is about 0.3 eV in water crystals, and even larger than 0.5 eV in ethylene crystals (even 
1.0 eV in (C2H 4 ) 10 ). These deviations can be decreased using larger fragmentation scheme: 
in REM' 8 , the typical deviations in water crystals decrease from 0.3 eV to about 0.15 eV, 
and in ethylene crystals, the typical deviations can reduce to about 0.1 eV. Obviously, such 
performances with BLYP and B3LYP are generally not satisfactory, since BLYP usually 
underestimate of non-local long-range electron-electron exchange interactions by the pure 
density functionals, and the B3LYP usually give wrong description of long-range interactions 
for this functional stem from modeling strong intermolecular interactions between solid 
macromolecular systems. This implies that the long-range corrections are very important 
for the REM-TDDFT calculations of large systems. 

TABLE V: Calculated S± excitation energies (in eV) by standard TDDFT and related 
excitation energy differences between REM-TDDFT results and them with different func- 
tionals in the ring molecular crystals 

BLYP/6-31+G* B3LYP/6-31+G* LC-BLYP/6-31+G* 

System 

TDDFT REM" 4 REM 5 TDDFT REM A REM B TDDFT REM A REM 5 

H2O ring 

(H 2 O)i 5.329 -0.266 -0.135 6.593 +0.113 -0.112 6.840 +0.076 -0.067 

(H 2 O) 20 5.355 +0.278 -0.155 6.626 +0.161 -0.143 6.837 +0.153 -0.054 

(H 2 O) 50 5.357 +0.336 -0.158 6.643 +0.155 -0.148 6.845 +0.144 -0.065 
C 2 H 4 ring 

(C 2 H 4 )i 5.329 -1.091 -0.338 6.297 -0.456 -0.135 7.108 -0.022 -0.071 

(C 2 H 4 ) 20 5.056 -0.847 -0.098 6.243 -0.445 -0.117 7.080 -0.031 -0.075 

(C 2 H 4 ) 50 5.038 -0.835 -0.090 6.235 -0.450 -0.120 7.074 -0.034 -0.080 

A REM-TDDFT with H 2 or C 2 H 4 as one monomer 

B REM-TDDFT with (H 2 0) 2 or (C 2 H 4 ) 2 as one monomer 

C. 2-D benzene crystal systems 

In this part, we attempt to apply the REM-TDDFT calculations to the 2-D benzene 
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crystal system. It is well known that crystals of acenes such as pentacene and tetracene own 
great potentials in the organic photovoltaic fields and consequently the accurate calculations 
of the electronic excited states of such aggregates are highly desired. The benzene crystal 
system can be seen as a simple model system of the acene crystals. Here we choose a number 
of benzene molecules in the (1,0,0) crystal face of the benzene crystal^ to do our test. As 
illustrated in Fig. [BJ each benzene column own its unique colour and four columns make 
up the 2-D benzene aggregates. One color square means one monomer and the dimers are 
automatically selected by the FMO subroutine in the GAMESS package. The LC-BLYP 
functional with 6-31G basis sets are used here. When doing the REM-TDDFT calculations, 
each monomer keeps one ground state (So) and one excited state (Si) and each dimer keeps 
one ground state (So) and two excited states (Si, S 2 ). 



The results of calculated Si excitation energies are listed in Table. [Vl] It could be found 
that the performance of REM-TDDFT exists in the 1-column situation: the difference be- 
tween the result of REM-TDDFT and TDDFT is only -0.003 eV. When the system tends to 
extend by adding the columns, the TDDFT results gradually converge to 5.672 eV. It means 
that the properties of the 2-D benzene system with such large sizes are already approaching 
the bulk ones of 2-D infinite benzene crystal. Here the results of REM-TDDFT are also con- 
verging very well to 5.660 eV, however, the difference for excitation energies of Si between 




OOOO 
# 



FIG. 8: The geometry of the 2-D benzene crystal. 
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REM-TDDFT and TDDFT increases to -0.012 eV, slightly larger than that in the 1-column 
case but still satisfactory. Such minor excitation energies errors with the magnitude from 
-0.003 eV to -0.012 eV for REM calculations of the 2-D benzene systems are comparable 
to those in the above 1-D examples, implying that REM has the potential to be applied to 
realistic molecular aggregates with complicated morphologies. 

TABLE VI: Calculated Si excitation energies (in eV) by 
REM-TDDFT and by standard TDDFT and the difference 
A (in eV) between the results of REM-TDDFT and TDDFT. 
The LC-BLYP/6-31G are used in the calculations. 



System 


1-Column 2-Columns 3-Columns 4-Columns 


REM-TDDFT 


5.682 


5.664 


5.660 


5.660 


TDDFT 


5.685 


5.671 


5.672 


5.672 


A 


-0.003 


-0.007 


-0.012 


-0.012 



D. Solvated systems 

Finally, let's turn our focus on the solutions. We choose two typical systems: one is the 
benzene + (H20) n system, the other is acetone + (H20) n system. These represent different 
solvation behaviours of non-polar and polar solutes dissolved in water. The geometries of 
the those two systems are chosen from the molecular dynamics (MD) trajectories in our 
other work—. Firstly, we performing the NVT MD simulations with simple point charge 
extended (SPC/E) potential for water and OPLS potential (optimized potentials for liquid 
simulations)^*^ for benzene and acetone, then select five uncorrelated snapshots for each 
set. From the selected snapshots, 12, 36 or 60 water molecules closest to benzene (acetone) 
as well as the central benzene (acetone) are taken to perform the REM-TDDFT calculations. 
Long-range corrected functional (LC-BLYP) with 6-31+G* basis functions are used here. 
There we treat one molecular unit as one fragment-monomer, then two molecular units as 
one fragment- dimer. In this test, we choose all of the benzene- water dimers, and the water- 
water dimers with interval less than 3. OA. Each monomer keeps the ground state So and 
the first singlet excited state Si, each dimer keeps So, Si and also the second singlet excited 
state S2. The results of both REM-TDDFT and standard TDDFT of these two solvated 
systems are listed in Table. IVHI and Table. IVIIIt separately. 
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The results of solvated benzene are listed in Table. I VIII It could be found that the 
REM-TDDFT can reproduce the full TDDFT values quite accurately. The average de- 
viation in the benzene + (H 2 0)i2 systems is only -0.006 eV, with mean square error of 
0.005 eV. The deviations will slightly increase with more water molecules: when adding 36 
H 2 Os, the average deviation increase to -0.014 eV, with mean square error of 0.011 eV; and 
when adding 60 H 2 Os, these two values turn to -0.027 eV and 0.021 eV, respectively. The 
deviations turn to large with the enlarged systems, for there are much more many-body 
interactions in those enlarged systems. Here we can also observe that the REM-TDDFT 
excitation energies are underestimated frequently relative to the full TDDFT values, the 
similar phenomenon is also observed in FM02-TDDFT calculations^. Introducing the 3- 
body interactions could improve the results , since only the two body interactions usually 
gives negative pair corrections^ 1 ^. 

The results of solvated acetone are listed in Table. IVIIIl It could be found that the REM- 
TDDFT can also reproduce the full TDDFT values quite well. The average deviations and 
the mean square errors (in brackets) are -0.020 eV (0.011 eV), -0.018 eV (0.073 eV) and - 
0.008 eV (0.066 eV), correspondingly. These deviations are larger than those in the solvated 
benzene, for the polar acetone is soluble in water, and there are stronger interactions between 
solute and solvent than those in solvated benzene systems. In solvated acetone systems, the 
electrons can be excited from acetone to acetone itself, and from acetone to the neighbor 
waters. While in the solvated benzene systems, the S\ excitations are mainly localized on the 
benzene molecule itself. In principle, one need to enlarge the fragment units or introducing 3- 
body (or higher) interactions to give a better description. Although lacking some interactions 
information, the wave functions of REM-TDDFT can also give qualitative correct pictures 
in these two solvated systems: the excitations are mainly from center benzene (or acetone) 
molecule (about 0.95-0.99, depends on the system); the contributions from waters are very 
small, and decrease with elongation of benzene-water (or acetone-water) distance. 

At last, we briefly introduce the timings for REM-TDDFT method. The time costs 
of both REM-TDDFT and standard TDDFT calculations for testing systems are listed in 
Table. IIX1 All calculations are implemented by the Sugon 12-core servers with Intel Xeon 
X5650@2.67GHz. It can be found from the table that the REM-TDDFT costs less time 
than standard TDDFT in this server, for the REM-TDDFT strategy has the approximate 
Ci x iV e 3 + n xc 2 x A^ 3 scaling when considering only the two-body interactions^. The 
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former is mainly from the calculation of the overlap matrix in the model space, and the 
latter is caused by the projections from target space to model space. There n is the number 
of dimers and iV e is the number of electrons of the whole system, while c\,Oi are constants 
affected by the preserved number of configurations in each state. Here the iV e 3 in mainly 



TABLE VII: Calculated S\ excitation energies (in eV) in solvated benzene systems by standard 
TDDFT and related excitation energy differences between REM-TDDFT results and them 



System 


benzene+(H20)i2 


benzene+ (H2 O ) 36 


benzene+(H20)6o 


REM TDDFT 


A* 


REM TDDFT 


A* 


REM TDDFT 


A* 


snapshot- 1 


5.349 


5.359 


-0.010 


5.322 


5.339 


-0.017 


5.317 


5.340 


-0.023 


snapshot-2 


5.309 


5.307 


0.002 


5.304 


5.323 


-0.019 


5.298 


5.355 


-0.057 


snapshot-3 


5.313 


5.317 


-0.004 


5.262 


5.291 


-0.029 


5.263 


5.296 


-0.033 


snapshot-4 


5.281 


5.286 


-0.005 


5.240 


5.250 


-0.010 


5.195 


.0 




snapshot-5 


5.268 


5.280 


-0.012 


5.272 


5.267 


0.005 


5.283 


5.280 


0.003 


A M 






_ _ 006 (0.005) 






-0.014(°- 011 ) 






_ a027 (0.021) 



* The difference is the result of REM-TDDFT minus that of standard TDDFT 
^ Accurate result is unavailable for the convergence problem in DFT 



TABLE VIII: Calculated S\ excitation energies (in eV) in solvated acetone systems by standard 
TDDFT and related excitation energy differences between REM-TDDFT results and them 



System 


acetone+ (H2 O) 12 


acetone-!- (H2 O ) 36 


acetone+ (H2 O )6o 


REM TDDFT 


A* 


REM TDDFT 


A* 


REM TDDFT 


A* 


snapshot- 1 


4.218 


4.229 


-0.011 


4.245 


4.274 


-0.029 


4.232 


4.290 


-0.058 


snapshot-2 


4.195 


4.210 


-0.015 


4.196 


4.227 


-0.031 


4.192 


4.229 


-0.037 


snapshot-3 


4.311 


4.345 


-0.034 


4.477 


4.354 


0.123 


4.488 


4.382 


0.106 


snapshot-4 


4.341 


4.374 


-0.033 


4.363 


4.430 


-0.067 


4.247 


_0 




snapshot-5 


4.368 


4.375 


-0.007 


4.340 


4.424 


-0.084 


4.366 


4.408 


-0.042 


A W 






-0.020(°- 011 ) 






_ aol8 (0.073) 






_ a008 (0.066) 



* The difference is the result of REM-TDDFT minus that of standard TDDFT 
^ Accurate result is unavailable for the convergence problem in DFT 
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from the lower triangular-upper triangular (LU) decomposition^, there the time-scale factor 
can at most up to N e , in practical applications it can be at most reduced to iV e Z d '°-2i. In 
fact, since REM-TDDFT using a disentanglement way^ 1 ^ to get the two body interactions, 
the various interactions extracting from dimers can be easily distribute to many servers, 
then the time costs can be even lower. 

TABLE IX: The approximate wall clock timing for 
REM-TDDFT and standard TDDFT calculations 
at LC-BLYP/6-31+G* level in the solvated acetone 
systems in Table. IVIII1 



System REM-TDDFT*' TDDFT 



acetone + (H 2 0)i2 


~0.5 


min 


~8 min 


acetone + (EbO^ 


~H 


min 


~200 min 


acetone + (H20)6o 


~58 


min 


~530 min 



* The timings of calculations of the monomers and 

dimers are not counted in. 
12-core server with Intel Xeon X5650 



V. SUMMARY AND CONCLUSION 

In this paper, we extend the ab initio REM method to TDDFT theory and use this ap- 
proach to calculate electronic excitation energies of various molecular aggregates. It is shown 
that this approach can not only gives a description of electronic excitation energies, but also 
provides a qualitative picture where the excitation locates. Since only the subsystems need 
to be solved in the whole aggregates, the computational costs are reduced remarkably than 
the TDDFT calculations while losing only little accuracy. Such achievements provide a 
new promising sub-system methodology for future quantitative studies of large complicated 
systems such as supramolecules, condensed phase matters. 

Test calculations for the one dimensional water molecule chains show that REM-TDDFT 
method is effective in reproducing the electronic excitation energies of low-lying excited 
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states: the typical deviation is only about 0.030 eV in Si or T\ states, and slight larger for 
higher excited states. The wave function analysis of REM-TDDFT also gives correct pic- 
tures of the excitation behavior in these systems. Furthermore, we test the REM-TDDFT 
with different basis sets and also various exchange-correlation functionals. We find that the 
larger basis sets will only slightly affect the final results, but the DFT functionals would 
significantly influence the stability and accuracy. Here the long-range corrected functionals 
with appropriate basis sets are recommended for dealing with large molecular aggregates. 
The trial test on the 2-D structure like benzene aggregates are also implemented and sat- 
isfactory excitation energy accuracies are also observed for them. At last, we turn to two 
types of aqueous systems to examine our REM-TDDFT's performances for the solutions. 
With LC-BLYP functional and 6-31 +G* basis sets, our REM-TDDFT method can repro- 
duce the standard TDDFT values quite well, for both of the aqueous systems with polar 
and non-polar solutes. 

The results of REM-TDDFT are acceptable in these molecular aggregate systems, how- 
ever, if one wants to pursue more accurate results the higher many-body interactions and 
ESP effects should be introduced. Progress along this direction is being made in our labo- 
ratory. 
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